{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2.2 Eigenvalue problems\n",
    "\n",
    "We solve the generalized eigenvalue problem\n",
    "\n",
    "$$\n",
    "A u = \\lambda M u,\n",
    "$$\n",
    "\n",
    "where $A$ comes from $\\int \\nabla u \\nabla v$, and $M$ from $\\int u v$, on the space $H_0^1$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import netgen.gui\n",
    "%gui tk\n",
    "from netgen.geom2d import unit_square\n",
    "from ngsolve import *\n",
    "from math import pi\n",
    "import scipy.linalg\n",
    "from scipy import random\n",
    "\n",
    "mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We setup a stiffness matrix $A$ and a mass matrix $M$, and declare a preconditioner for $A$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "fes = H1(mesh, order=4, dirichlet=[1,2,3,4])\n",
    "u = fes.TrialFunction()\n",
    "v = fes.TestFunction()\n",
    "\n",
    "a = BilinearForm(fes)\n",
    "a += SymbolicBFI(grad(u)*grad(v))\n",
    "pre = Preconditioner(a, \"multigrid\")\n",
    "\n",
    "m = BilinearForm(fes)\n",
    "m += SymbolicBFI(u*v)\n",
    "\n",
    "a.Assemble()\n",
    "m.Assemble()\n",
    "\n",
    "u = GridFunction(fes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The inverse iteration is\n",
    "\n",
    "$$\n",
    "u_{n+1} = A^{-1} M u_n,\n",
    "$$\n",
    "\n",
    "where the Rayleigh quotient\n",
    "$$\n",
    "\\rho_n = \\frac{\\left \\langle A u_n, u_n\\right \\rangle}{\\left \n",
    "\\langle M u_n, u_n\\right \\rangle}\n",
    "$$\n",
    "\n",
    "converges to the smallest eigenvalue, with rate of convergence $\\lambda_1 / \\lambda_2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The preconditioned inverse iteration (PINVIT), see [Knyazef+Neymeyr], replaces $A^{-1}$ by an approximate inverse $C^{-1}$:\n",
    "\n",
    "$$\n",
    "\\rho_n = \\frac{\\left \\langle A u_n, u_n\\right \\rangle}{\\left \\langle M u_n, u_n\\right \\rangle} \\\\\n",
    "w_n = C^{-1} (A u_n - \\rho M u_n) \\\\\n",
    "u_{n+1} = u_n + \\alpha w_n\n",
    "$$\n",
    "\n",
    "The optimal step-size $\\alpha$ is found by minimizing the Rayleigh-quotient on a two-dimensional space:\n",
    "\n",
    "$$\n",
    "u_{n+1} = \\operatorname{arg} \\min_{v \\in \\{ u_n, w_n\\}} \\frac{\\left \\langle A v, v\\right \\rangle}{\\left \\langle M v, v\\right \\rangle} \n",
    "$$\n",
    "\n",
    "This minimization problem can be solved by a small eigenvalue problem\n",
    "\n",
    "$$\n",
    "a y = \\lambda m y\n",
    "$$\n",
    "with matrices\n",
    "$$\n",
    "a = \n",
    "\\left( \\begin{array}{cc}\n",
    "\\left \\langle  A u_n, u_n \\right \\rangle &\n",
    "\\left \\langle  A u_n, w_n \\right \\rangle \\\\\n",
    "\\left \\langle  A w_n, u_n \\right \\rangle &\n",
    "\\left \\langle  A w_n, w_n \\right \\rangle \n",
    "\\end{array} \\right) \\quad\n",
    "m = \n",
    "\\left( \\begin{array}{cc}\n",
    "\\left \\langle  M u_n, u_n \\right \\rangle &\n",
    "\\left \\langle  M u_n, w_n \\right \\rangle \\\\\n",
    "\\left \\langle  M w_n, u_n \\right \\rangle &\n",
    "\\left \\langle  M w_n, w_n \\right \\rangle \n",
    "\\end{array} \\right)\n",
    "$$\n",
    "\n",
    "Then, the new iterate is\n",
    "\n",
    "$$\n",
    "u_{n+1} = y_1 u_n + y_2 w_n\n",
    "$$\n",
    "\n",
    "where $y$ is the eigenvector corresponding to the smaller eigenvalue."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementation in NGSolve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create some help vectors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = u.vec.CreateVector()\n",
    "w = u.vec.CreateVector()\n",
    "Mu = u.vec.CreateVector()\n",
    "Au = u.vec.CreateVector()\n",
    "Mw = u.vec.CreateVector()\n",
    "Aw = u.vec.CreateVector()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "random initial condition, zero on the boundary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "r.FV().NumPy()[:] = random.rand(fes.ndof)\n",
    "u.vec.data = Projector(fes.FreeDofs(), True) * r\n",
    "\n",
    "# from random import random\n",
    "# freedofs = fes.FreeDofs()\n",
    "# for i in range(len(u.vec)):\n",
    "#    u.vec[i] = random() if freedofs[i] else 0  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "26.310136660400275\n",
      "2.0782284240636826\n",
      "2.0016788490372455\n",
      "2.0002002507016474\n",
      "2.0000413285947802\n",
      "2.0000105080928825\n",
      "2.0000027530788227\n",
      "2.0000007313313577\n",
      "2.0000001959000486\n",
      "2.0000000528126267\n",
      "2.0000000143203\n",
      "2.000000003918371\n",
      "2.0000000010972063\n",
      "2.000000000329983\n",
      "2.0000000001208567\n",
      "2.0000000000637526\n",
      "2.000000000048131\n",
      "2.000000000043857\n",
      "2.0000000000426827\n",
      "2.0000000000423617\n"
     ]
    }
   ],
   "source": [
    "for i in range(20):\n",
    "    Au.data = a.mat * u.vec\n",
    "    Mu.data = m.mat * u.vec\n",
    "    auu = InnerProduct(Au, u.vec)\n",
    "    muu = InnerProduct(Mu, u.vec)\n",
    "    # Rayleigh quotiont\n",
    "    lam = auu/muu\n",
    "    print (lam / (pi**2))\n",
    "    # residual\n",
    "    r.data = Au - lam * Mu\n",
    "    w.data = pre.mat * r.data\n",
    "    w.data = 1/Norm(w) * w\n",
    "    Aw.data = a.mat * w\n",
    "    Mw.data = m.mat * w\n",
    "\n",
    "    # setup and solve 2x2 small eigenvalue problem\n",
    "    asmall = Matrix(2,2)\n",
    "    asmall[0,0] = auu\n",
    "    asmall[0,1] = asmall[1,0] = InnerProduct(Au, w)\n",
    "    asmall[1,1] = InnerProduct(Aw, w)\n",
    "    msmall = Matrix(2,2)\n",
    "    msmall[0,0] = muu\n",
    "    msmall[0,1] = msmall[1,0] = InnerProduct(Mu, w)\n",
    "    msmall[1,1] = InnerProduct(Mw, w)\n",
    "    # print (\"asmall =\", asmall, \", msmall = \", msmall)\n",
    "    \n",
    "    \n",
    "    eval,evec = scipy.linalg.eigh(a=asmall, b=msmall)\n",
    "    # print (eval, evec)\n",
    "    u.vec.data = float(evec[0,0]) * u.vec + float(evec[1,0]) * w\n",
    "    \n",
    "Draw (u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simultaneous iteration for several eigenvalues"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Declare GridFunction with multiple components to store several eigenvectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "num = 5\n",
    "u = GridFunction(fes, multidim=num)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create list of help vectors, random initialization of u vectors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = u.vec.CreateVector()\n",
    "Av = u.vec.CreateVector()\n",
    "Mv = u.vec.CreateVector()\n",
    "\n",
    "vecs = []\n",
    "for i in range(2*num):\n",
    "    vecs.append (u.vec.CreateVector())\n",
    "\n",
    "for v in u.vecs:\n",
    "    r.FV().NumPy()[:] = random.rand(fes.ndof)\n",
    "    v.data = Projector(fes.FreeDofs(), True) * r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute num residuals, and solve small eigenvalue problem on $2 num$-dimensional space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 : [11.6130455768801, 61.066110947635224, 82.49711365113266, 90.41384577806575, 96.98583176269047]\n",
      "1 : [2.037296014502914, 5.994129376899316, 9.212270101228079, 16.752145334653957, 19.356324864358374]\n",
      "2 : [2.0014100615066996, 5.092284538737162, 5.238911096053245, 10.123988712254931, 10.879020322598242]\n",
      "3 : [2.0002554163633968, 5.012268562196962, 5.016722138616416, 9.213915664378316, 10.256909163977726]\n",
      "4 : [2.000057169527135, 5.001025982274766, 5.0025582296087965, 8.719831456835285, 10.0918178567872]\n",
      "5 : [2.000013593722611, 5.000120805946273, 5.000384594365776, 8.377992268553804, 10.035379830359405]\n",
      "6 : [2.000003314919647, 5.000020247177596, 5.000055876284345, 8.184287429684003, 10.01410422992045]\n",
      "7 : [2.000000821205276, 5.000004016939583, 5.00000849353781, 8.08609184849783, 10.005671478126507]\n",
      "8 : [2.0000002058391035, 5.0000008179428255, 5.000001478437589, 8.039447631058222, 10.002288011579898]\n",
      "9 : [2.0000000520461136, 5.000000168593007, 5.000000304198153, 8.017928866313158, 10.000924591937792]\n",
      "10 : [2.000000013279318, 5.000000039218533, 5.000000071766628, 8.008121745129989, 10.000374075767763]\n",
      "11 : [2.0000000034302508, 5.000000012511654, 5.000000020495333, 8.003675600991937, 10.000151525029434]\n",
      "12 : [2.0000000009138295, 5.000000006632997, 5.0000000086298915, 8.00166315351873, 10.00006147486362]\n",
      "13 : [2.000000000267371, 5.000000005221903, 5.000000005864453, 8.00075273753271, 10.000025013254238]\n",
      "14 : [2.0000000001006097, 5.000000004822398, 5.000000005260717, 8.00034080330491, 10.000010241080249]\n",
      "15 : [2.000000000057419, 5.00000000470997, 5.000000005129163, 8.000154364958556, 10.000004253440505]\n",
      "16 : [2.0000000000461995, 5.000000004681857, 5.000000005098049, 8.00006995469134, 10.000001825435998]\n",
      "17 : [2.000000000043275, 5.000000004674421, 5.000000005090545, 8.000031724949128, 10.000000840538995]\n",
      "18 : [2.0000000000425104, 5.000000004672643, 5.000000005088353, 8.000014405452335, 10.000000440900093]\n",
      "19 : [2.00000000004231, 5.000000004672329, 5.000000005087833, 8.000006556850161, 10.000000278698185]\n"
     ]
    }
   ],
   "source": [
    "asmall = Matrix(2*num, 2*num)\n",
    "msmall = Matrix(2*num, 2*num)\n",
    "lams = num * [1]\n",
    "\n",
    "for i in range(20):\n",
    "    \n",
    "    for j in range(num):\n",
    "        vecs[j].data = u.vecs[j]\n",
    "        r.data = a.mat * vecs[j] - lams[j] * m.mat * vecs[j]\n",
    "        vecs[num+j].data = pre.mat * r\n",
    "\n",
    "    for j in range(2*num):\n",
    "        Av.data = a.mat * vecs[j]\n",
    "        Mv.data = m.mat * vecs[j]\n",
    "        for k in range(2*num):\n",
    "            asmall[j,k] = InnerProduct(Av, vecs[k])\n",
    "            msmall[j,k] = InnerProduct(Mv, vecs[k])\n",
    "\n",
    "    ev,evec = scipy.linalg.eigh(a=asmall, b=msmall)\n",
    "    lams[:] = ev[0:num]\n",
    "    print (i, \":\", [lam/pi**2 for lam in lams])\n",
    "    \n",
    "    for j in range(num):\n",
    "        u.vecs[j][:] = 0.0\n",
    "        for k in range(2*num):\n",
    "            u.vecs[j].data += float(evec[k,j]) * vecs[k]\n",
    "            \n",
    "Draw (u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The *multidim-component* select in the *Visualization* dialog box allows to switch between the components of the multidim-GridFunction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
